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Abstract. The Bayesian approach has proved to be a coherent approach to handle iU 
posed Inverse problems. However, the Bayesian calculations need either an optimization 
or an integral calculation. The maximum a posteriori (MAP) estimation requires the 
minimization of a compound criterion which, in general, has two parts: a data fitting 
part and a prior part. 

In many situations the criterion to be minimized becomes multimodal. The cost of 
the Simulated Annealing (SA) based techniques is in general huge for inverse problems. 
Recently a deterministic optimization technique, based on Graduated Non Convexity 
(GNC), have been proposed to overcome this difficulty. 

The objective of this paper is to show two specific implementations of this technique for 
the following situations: 

- Linear inverse problems where the solution is modeled as a piecewise continuous func- 
tion. The non convexity of the criterion is then due to the special choice of the prior; 

- A nonlinear inverse problem which arises in inverse scattering where the non convexity 
of the criterion is due to the likelihood part. 

Key words: Inverse problems, Regularization, Bayesian calculation. Global optimiza- 
tion. Graduated Non Convexity 

1. Introduction 

We consider tlie case of general inverse problems: 



where x is the vector of unknown variables, y is the data, A is a linear or non 
linear operator and b represents the errors which are assumed, hereafter, additive, 
zero-mean, white and Gaussian. 

The Bayesian approach has proved to be a coherent approach to handle these 
problems. However, the Bayesian calculations need either an optimization or an 
integral calculation. The maximum a posteriori (MAP) estimation needs the min- 
imization of a compound criterion: 



y = A{x) + b, 





(2) 



which, in general, has two parts: a data fitting part: 

Q{x) = -logp{y\x) = \\y - Aix)f 



(3) 
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and a prior part: 

This last expression is due to general Markov modeling where r and s are two site 
indexes and Afr means neighbor sites of r. 

In many situations the criterion J{x) becomes multimodal. We consider here 
two cases: The first is the case of general linear inverse problems with markovian 
priors with non convex energies, and the second is, the case of non- linear inverse 
problem. 

In both cases we need a global optimization technique to determine the solution. 
In the first case the non convexity is due to the second termed and in the second 
case the non convexity is more due to the first termed. 

The cost of the Simulated Annealing (SA) based techniques is mainly depen- 
dent to the neighborhood size of the posterior marginal probability distribution 
p{xj\x^ y) which is directly related to the neighborhood size of the prior marginal 
probability distribution p{xj\x) and the support of the operator A. When A is 
a local operator with very small support, i.e.; when the data element yi depends 
only to a few number of the unknown variables Xj, then SA can be implemented 
efhciently 1). But, unfortunately, this is not the case of many inverse prob- 

lems, where the support of the operator is not small. The cost of SA is then, in 
general, huge for these problems. 

Recently a deterministic relaxation algorithm, inspired by the Graduated Non 
Convexity (GNC) principle, has been proposed by Blake and Zisserman in |5| for 
the optimization of the multimodal MAP criteria. They have shown its efficiency in 
practical applications for noise cancellation and segmentation. This algorithm has 
been extended to the general linear ill-posed inverse problem by Nikolva et aL in Q . 

The object of this presentation is to show two specific implementations of this 
technique for two specific cases of the two aforementioned situations, i.e.; 

- The linear inverse problems where the solution is modeled as a piecewise contin- 
uous function using a compound markov modeling (for example the intensity and 
the line process in image reconstruction), where the non convexity of the criterion 
is due to the markovian priors; and 

- A special non-linear inverse problem which arises in inverse scattering and diffrac- 
tion tomography imaging applications where the non convexity of the criterion is 
due to the likelihood part. 

The paper is organized as follows: The next section presents the main idea of 
the GNC principle. Sections 3 and 4 will consider the two aforementioned specific 
cases and, finally, some simulation results will illustrate the performances of the 
proposed method in two special applications. 
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2. Graduated Non Convexity scheme 

The principle of this algorithm is very simple. It consists of approximating the 
non convex criterion J{x) with a sequence of continuously derivable criteria Jc^ (x) 
such that: 

• the first one Jc,, i^) be convex; 

• the final one (the limit) Jcki^) converges to J{x): 

lim Jck{x) — J{x), Wx, 

where Ck > are increasing relaxation parameters; and then 

• for each fc, a relaxed solution is calculated by minimizing locally, initialized by 
the previous solution, as follows: 

2cfc=arg mm Udx)} (5) 
xev{x,^^^) 

Fig. ^ illustrates such a scheme. The hope is that the sequence Xc^, converges to 
the global minimizer of J{x). Note that there is no theoretical ground for this 
hope, however, in many practical applications it seems to be realist. 




Figure 1: GNC scheme. 
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3. Lineeir inverse problems with a piecewise Gaussian prior 

Let first consider a very simple noise filtering problem: 

y = x + b, (6) 

where we know that x represents the samples of a piecewise continuous function 
x{t). Blake and Zisserman in [Q proposed to estimate x by searching the global 
minimum of the following criterion: 

J{x)^\\y~xr + n{x) (7) 

with 

j 

and 

^^^> '{ a = {\Tf if \t\ > T 
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Figure 2: Truncated quadratic function. 

Note that, for (pit) — the criterion J{x) can be considered either as the 
MAP criterion with Gaussian prior or as the Tikhonov regularization one. The 
choice (H) for (j){t) is done to preserve the discontinuities in x{t). With this choice, 
obviously, J{x) is multimodal. The GNC idea was then to construct: 

J,{x) = \\y-xW' + n,{x) (10) 

with 

n,(x)^Y.'i^,{t,), (11) 

and 

r {\tf if |i| < q, 

0c(t) = < a-\c{\t\^r,f if qc>\t\>r, 
y a = (AT)2 if \t\ > 



gc = r(l + 2AVc)-i/2 
rc = T(l + 2AVc)i/2 



(12) 
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and find a cq such that Jc for c < co be convex and then for a given sequence of 
relaxations parameters {cq, ci, . . . , Cfc} do: 

2c, =arg min {Jcu{x)} (13) 
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Figure 3: Truncated quadratic function and it's relaxation scheme. 
Blake and Zisserman in Q showed the existence of cq such that Jcoi^) be 
convex and global convergence of this algorithm. However, for the case of inverse 
problems y = Ax + b where A is singular or ill-conditioned, the existence of Cg is 
no more insured. Nikolova et al. extended this work by proposing a doubly 

relaxed criterion: 

Jc{x) = \\y - x\\^ + ilaAx) (14) 

with 

f^a,c(a;) =^[0c(tj)+ai'] (15) 
j 

and the following double relaxation scheme: 

for fixed c = cq and for a — oq, . . . , do: 

2a, =arg niin {Ja,,coix)} (16) 
xev{x^^_^) 

for fixed a = and for c — cq, . . . ,oo do: 

Xc^ = arg mm {Jc^Ax)} (17) 

xev{x,^_^) 

The initial convexity of the criterion is insured for > Many details, discus- 
sions and more extensions, specially for 2D case, are given in |^, ^. 
3.1. Link with compound Markov models 

Compound Markov modeling in image processing became popular after the 
works of Geman & Geman Q and Besag 1 10| ll| . To see the link with these model 



briefly, let consider a 1-D case. When x{t) is assumed piecewise continuous or 
piecewise Gaussian, it can be modeled with a set of coupled variables [x, I), where 
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Z is a vector of binary- valued variables and x is Gaussian. Now, consider the MAP 
estimate of {x,l): 

(xj) = argmin {J{x,l) = - logp(a;, Z|y)} (18) 
{x,l) 

where 

- ^ogp{x, l\y) = - logp{y\x) - logp{x\l) - logp{l) + cte, (19) 
with the following prior laws: 

-logp(y|a;) = ||y-Aa?f 

-logp(a;|Z) = ^0(t^.)(i_;^.), tj=xj-xj^i, cj)[t) = {\tf 







-logp(Z) = 

j 

we obtain: 

{xA) = arg min {J{x. I)} (20) 
ix,l} 

with 

Jix, I) = \\y - Axf + [>^Hxj - x,^i)^{l - I,) + alj] . (21) 

j 

Note that the line variables Ij are assumed non-interacting (mutually independent). 
With this hypothesis, it is easy to show that the solution {x,l) obtained by ( pO| ) 
is equivalent to the solution obtained by 

S = argmin {J(a;)} with J{x) = \\y - AxW"^ + 'Y^(j>{xj — xj^i)'^ (22) 

j 

with (/)(t) defined by and 

h = { I fj-l^^'l'^l (23) 
^ \ if \xj - Xj^i \ <T ^ ^ 



4. A non linear inverse problem 

To show how the GNC principle can be used for nonlinear inverse problems we 
consider the case of inverse scattering and more specifically the diffraction tomog- 
raphy. To be short in presentation of the application, we give here an abstract 
presentation of the problem. For more details on derivation of the inverse scat- 
tering and diffraction tomography application see |l^, |l3|, |l^. To summarize, 
considering the geometry of Fig. ^ we have the following relations: 

y{n) = JJ Gm{r,,r')(l>ir')x{r')dr', r, e S, (24) 
0(r) = Mr)+ jj Go{r,r')(l){r')x{r')Ar', r e D. (25) 
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The discretized version of these two equations can be written with the following 
compact notations: 

y = G„^Xcj) + b, (26) 
</. = ct,o + GoXcf>, (27) 

where 4>q is the incident field, G„i, Go are matrices related to the Green functions, 
is a diagonal matrix with the components of the vector x (a, n length vector) as 
its diagonal elements and y, (p are respectively m and n length vectors representing 
the measured data (scattered field) and the total field on the object. Note that n 
may be greater than m. 




Figure 4: Diffraction tomography geometry configuration 

These two equations can be combined to obtain a symbolic explicit relation 
between the data y and the unknowns x: 



y = G,„X (/ - GoX) ^(po + b = Aix) + b, 



(28) 



where the considered matrix is assumed to be invertible. Now, the inverse problem 
we are faced is to find x given y. Note that, (pQ, Go and G„i are known and the 
relation between x and y is non linear. In fact, given (j), y is linear in x (p6|), but 
(p depends on x through the second equation (p7|). 

Here also, using the Bayesian approach, the MAP estimate is defined as the 
minimize! of 

J{x) ^\\y - Aix)f + n{x) (29) 

and, even when fl{x) is chosen to be convex, J{x) may not due to the fact that 
\\y — A{x)\\'^ is no more quadratic in x. Carfantan et al. in ||T^, |T^, I4| proposed 
to use the GNC idea in this case. 

To introduce the GNC technique, they considered the following relaxation se- 
quence: 

Ao, {x) = GmXil - ckGoX)-^(t>o. (30) 

with Co = 0, and lim^-^oo Cfe = I. 

Note that the first term (co = 0) corresponds to a linearized model for the 
problem named the Born approximation which consists in neglecting partially the 
diffraction effects. This results to the following convex criterion: 

Jo{x) = ||y-G™X0of + 17(a;). 
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Note also that for = 1 the criterion Ji(x) = J{x). The main practical problem 
is then the choice of sequences {ck = 0, • • • , 1} which is done by experiment. 



5. Some simulation results and applications 

5.1. 1-D NOISE FILTERING 



Fig. 5.1. shows an example of results obtained in noise filtering. In this figure 
we see the original signal, noisy data, restoration using a Gaussian model {4>{t) 
quadratic) and restoration obtained by GNC when (f){t) is chosen to be truncated 
quadratic. 



Original & GNC reconstruction 
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Figure 5: Noise filtering: 

a) Original and data, b) Original, Gaussian restoration and GNC restoration 

5.2. 1-D SIGNAL DECONVOLUTION 



Fig. p.2j shows an example of results in signal deconvolution. In this figure 
we see the original signal, noisy data, restoration using a Gaussian model and 
restoration obtained by GNC. 



5.3. Image restoration 



Fig. 5.3. shows an example of results in Image restoration. In this figure a) 
is the original image, b) is the blurred and noisy data, c) is the restoration using 
a Gaussian model and d) is the restoration obtained by GNC with truncated 
quadratic regularization. 
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Original & GNC reconstruction 




Figure 6: Deconvolution: 

a) Original and data, b) Original, Gaussian restoration and GNC restoration 
5.4. Image reconstruction in X-ray tomography 



Fig. 5.4. shows results obtained in X-ray tomography image reconstruction. In 
this figure a) shows the original image, b) shows the projections (data), c) shows 
the backprojection reconstruction, d) shows a reconstruction using a Gaussian 
prior, e) shows a reconstruction using a Gamma prior, and f) shows a reconstruc- 
tion using GNC with truncated quadratic regularization. 



5.5. Inverse scattering and diffraction tomography 



Fig. p.5.| shows an example of results in non linear diffraction tomography image 
reconstruction. In this figure a) is the original image, b) is the measured scattered 
field data, d) is a reconstruction using the linear Born approximation, and e) is a 
reconstruction using GNC. 



6. Conclusions 

The Bayesian maximum a posteriori (MAP) estimates requires the minimization 
of a compound criterion which, in general, has two parts: a data fitting part and 
a prior part. In many situations in inverse problems the criterion to be minimized 
is multimodal. The cost of the Simulated Annealing (SA) based techniques is in 
general huge for these problems. 

We reported here recently proposed new techniques, based on Graduated Non 
Convexity (GNC), to overcome this difficulty and showed two specific implemen- 
tations of this technique for: 

- The linear inverse problems such as: noise filtering, deconvolution, image restora- 
tion and tomographic image reconstruction, where the solution is modeled as a 
piecewise continuous function and where, the non convexity of the criterion is due 
to the special choice of the prior; and 

- A nonlinear inverse inverse scattering and diffraction tomography where the non 
convexity of the criterion is due to the likelihood part. 



10 



A. Mohammad-Djafari et al. 



References 

1. S. Geman and D. Geman, "Stochastic relaxation, Gibbs distributions, and the 
Bayesian restoration of images," IEEE Transactions on Pattern Analysis and Ma- 
chine Intelligence, vol. PAMI-6, pp. 721-741, Nov. 1984. 

2. L. Younes, "Estimation and annealing for Gibbsian fields," Annales de I'tnstitut 
Henri Poincare, vol. 24, pp. 269-294, Feb. 1988. 

3. F. Jeng and J. Woods, "Simulated annealing in compound Gaussian random fields," 
IEEE Transactions on Information Theory, vol. IT-36, pp. 94-107, Jan. 1990. 

4. A. Blake and A. Zisserman, Visual reconstruction. Cambridge: The MIT Press, 
1987. 

5. A. Blake, "Comparison of the efiiciency of deterministic and stochastic algorithms 
for visual reconstruction," IEEE Transactions on Pattern Analysis and Machine 
Intelligence, vol. PAMI-11, pp. 2-12, January 1989. 

6. M. Nikolova, A. Mohammad-Djafari, and J. Idler, "Inversion of large-support ill- 
conditioned linear operators using a Markov model with a line process," in ICASSP, 
vol. V, (Adelaide, Austraha), pp. 357-360, 1994. 

7. M. Nikolova and A. Mohammad-Djafari, "Discontinuity reconstruction from linear 
attenuating operators using the weak-string model," in Proceedings of European Sig- 
nal Processing. Conf, vol. 2, pp. 1062 1066, 1994. 

8. M. Nikolova, Inversion markovienne de problemes lineaires mal poses, application a 
I'imagerie tomographique. PhD thesis, Universite de Paris-Sud, Orsay, Feb. 1995. 

9. M. Nikolova, J. Idler, and A. Mohammad-Djafari, "Inversion of large-support ill- 
posed linear operators using a piecewise Gaussian mrf," tech. rep., GPI-LSS, sub- 
mitted to IEEE Transactions on Image Processing, Gif-sur-Yvcttc, France, 1995. 

10. J. E. Besag, "Spatial interaction and the statistical analysis of lattice systems (with 
discussion)," Journal of the Royal Statistical Society B, vol. 36, no. 2, pp. 192-236, 
1974. 

11. J. E. Bcsag, "On the statistical analysis of dirty pictures (with discussion)," Journal 
of the Royal Statistical Society B, vol. 48, no. 3, pp. 259-302, 1986. 

12. H. Carfantan and A. Mohammad-Djafari, "A Bayesian approach for nonlinear in- 
verse scattering tomographic imaging," in ICASSP, vol. IV, (Detroit, U.S.A.), 
pp. 2311-2314, May 1995. HC. 

13. H. Carfantan and A. Mohammad-Djafari, "Approche bayesiennc ct algorithmc mul- 
tiresolution pour un problcmo inverse non lineaire en tomographic dc diffraction," in 
Actes du 15^ Colloque GRETSI, vol. 2, (Juan-les -pins, France), pp. 849-852, Sept. 
1995. 

14. H. Carfantan and A. Mohammad-Djafari, "Beyond the Born approximation in in- 
verse scattering with a Bayesian approach," in 2nd International Conference on 
Inverse Problems in Engineering, (Le Croisic, France), June 1996. 



NEW ADVANCES IN BAYESIAN CALCULATION. 



11 



30 40 

a 



10 20 30 40 50 60 

b 



10 20 30 40 50 60 



10 20 30 40 50 60 



Figure 7: Image restoration: 

a) original, b) data, c) Gaussian restoration, d) GNC restoration 
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Figure 9: Diffraction tomography image reconstruction: 

a) original, b) linear Born approximation reconstruction, c) GNC based recon- 
struction 



